#include "VGQSeries.h"
#include <math.h>
#include <gsl\gsl_sf_gamma.h>
#include <boost/math/special_functions/factorials.hpp>
#include <boost/math/special_functions/bessel.hpp>
#include <boost/math/special_functions/binomial.hpp>

VGQSeries::VGQSeries(double mu, double lambda, double alpha, double beta) 
	: _mu(mu), _lambda(lambda), _alpha(alpha), _beta(beta), 
		_gamma(sqrt(alpha*alpha - beta*beta)), _vq(_coefficients),
		_pvgDistribution(new VarianceGammaDistribution(0.0, _lambda, _alpha, _beta))
{
	if(!_pvgDistribution->validParameters(mu, lambda, alpha, beta))
		throw std::exception("VG parameters invalid");

	initialize();
}

VGQSeries::VGQSeries( std::shared_ptr<VarianceGammaDistribution > pvgDistribution ) : _pvgDistribution(pvgDistribution), _vq(_coefficients)
{
	_mu = _pvgDistribution->mu();
	_lambda = _pvgDistribution->lambda();
	_alpha = _pvgDistribution->alpha();
	_beta = _pvgDistribution->beta();
	_gamma = _pvgDistribution->gamma();

	initialize();
}

void VGQSeries::initialize()
{
	_N0 = 1/_pvgDistribution->normalizingConstant();

	_zeroQuantileLocation = _pvgDistribution->zeroQuantileLocation();
}

VGQSeries::~VGQSeries(void)
{
}

double VGQSeries::A( unsigned int n, double y )
{
	return exp(-_beta * ( y + _x0 )) * pow(-_beta, (int)n);
}

double VGQSeries::a( unsigned int n )
{
	using namespace boost::math;
	size_t size = _va.size();
	if(n < size)
		return _va[n];
	
	_va.resize(n + 1, 0.0);

	for(unsigned int i = size; i < _va.size(); i++)
	{
		_va[i] = A(i, 0) / factorial<double>(i);
	}

	return _va[n];
}

double VGQSeries::B( unsigned int n, double y )
{
	double poch = (n==0) ? 1 / (_lambda - 0.5) : gsl_sf_poch(0.5 + _lambda, n - 1);

	if(y >= -_x0)
		return ( std::pow(-1.0, int(n)) * std::pow(_x0 + y, 0.5 - n - _lambda)*(-1 + 2*_lambda)* poch )/2;
	
	return (pow(-_x0 - y,0.5 - n - _lambda)*(-1 + 2*_lambda) * poch)/2;
}

double VGQSeries::b( unsigned int n )
{
	using namespace boost::math;

	size_t size = _vb.size();
	if(n < size)
		return _vb[n];

	_vb.resize(n + 1, 0.0);

	for(unsigned int i = size; i < _vb.size(); i++)
	{
		_vb[i] = B(i, 0) / factorial<double>(i);
	}

	return _vb[n];
}

double VGQSeries::C( unsigned int n, double y )
{
	using namespace boost::math;
		
	int sign = (y >= -_x0) ? 1 : -1;
	double sum = 0.0;
	for(unsigned int k = 0; k <= n; k++)
	{
		 sum += binomial_coefficient<double>(n, k) * cyl_bessel_k(_lambda - 0.5 - (2*k) + n, sign * _alpha * (y + _x0));	
	}

	return std::pow((double)sign, (int)n) * std::pow(-_alpha/2.0,(int)n) * sum;
}

double VGQSeries::c( unsigned int n )
{
	using namespace boost::math;

	size_t size = _vc.size();
	if(n < size)
		return _vc[n];

	_vc.resize(n + 1, 0.0);

	for(unsigned int i = size; i < _vc.size(); i++)
	{
		_vc[i] = C(i, 0) / factorial<double>(i);
	}

	return _vc[n];
}

double VGQSeries::d( unsigned int n )
{
	if(n < _vd.size())
		return _vd[n];

	if(n == 0){
		_vd.push_back(1/c(0));
		return _vd[0];
	}

	double sum = 0.0;
	for(unsigned int k = 1; k <= n; k++)
	{
		sum += c(k)*d(n-k);
	}

	_vd.resize(n+1, 0.0);
	_vd[n] = (-1/c(0))*sum;

	return _vd[n];
}

double VGQSeries::g( unsigned int n )
{
	size_t size = _vg.size();
	if(n < size)
		return _vg[n];

	_vg.resize(n + 1, 0.0);

	for(size_t i = size; i < _vg.size(); i++)
	{
		double sum = 0.0;

		for(size_t k = 0; k <= i; k++)
		{
			for(size_t j = 0; j <= k; j++)
			{
				sum += ( d(i-k) * b(k-j) * a(j));  
			}
		}

		_vg[i] = _N0*sum;
	}

	return _vg[n];
}

double VGQSeries::h( unsigned int n )
{
	size_t size = _vh.size();
	if(n < size)
		return _vh[n];

	if(n == 0){
		_vh.push_back(g(0));
		return _vh[0];
	}

	double sum = 0.0;
	for(size_t k = 0; k <= n; k++)
	{
		sum += ( g(k) * q(n, k) );			
	}
	
	_vh.resize(n + 1, 0.0);
	_vh[n] = sum;

	return _vh[n];
}

double VGQSeries::q( unsigned int n, unsigned int k )
{
	auto it = _mq.find(k);

	if(it == _mq.end())
		it = _mq.insert(pair<size_t, vector<double> >(k, vector<double>())).first;
		
	auto& vqk = it->second;
	if(n < vqk.size())
		return vqk[n];

	if(n < k){
		return 0.0;
	}

	if(n == k){
		vqk.resize(n+1,0.0);
		vqk[n] = pow(q(1), (int)k);
		return vqk[n];
	}

	if(k == 0)
		if(n==0)
			return 1.0;
		else
			return 0.0;

	if(k==1)
		return q(n);

	double sum = 0.0;
	for(size_t j=1; j <= (n-k); j++)
	{
		double term2 = ((double(k+1)/double(n-k))*double(j))-1.0;
		if(term2 == 0)
			continue;
		
		double term1 = q(n-j,k);
		if(term1 == 0)
			continue;

		sum += term1*term2*(q(j+1)/q(1));
	}

	vqk.resize(n+1,0.0);
	vqk[n] = sum;
		
	return vqk[n];
}

double VGQSeries::q( unsigned int n )
{
	size_t size = _vq.size();

	if(n < size)
		return _vq[n];

	if(n == 0){
		_vq.push_back(0.0);
		return _vq[n];
	}

	double val = (1/double(n))*h(n-1);

	_vq.resize(n + 1, 0.0);
	_vq[n] = val;

	return _vq[n];
}

void VGQSeries::SetInitialCondition( double u )
{
	clearCache();

	_u0 = u;
	_x0 = _pvgDistribution->quantile(_u0);
}

void VGQSeries::clearCache()
{
	_va.clear();
	_vb.clear();
	_vc.clear();
	_vd.clear();
	_vg.clear();
	_vh.clear();
	_vq.clear();
	_mq.clear();
}



